Thermal diffusion by Brownian motion induced fluid stress 
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The Ludwig-Soret eflect, the migration of a species due to a temperature gradient, has been 
extensively studied without a complete picture of its cause emerging. Here we investigate the 
dynamics of DNA and spherical particles subjected to a thermal gradient using a combination of 
Brownian dynamics and the lattice Boltzmann method. We observe that the DNA molecules will 
migrate to colder regions of the channel, an observation also made in the experiments of Duhr, et 
al[3|. In fact, the thermal diffusion coefficient found agrees quantitatively with the experimental 
value. We also observe that the thermal diffusion coefficient decreases as the radius of the studied 
spherical particles increases. Furthermore, we observe that the thermal fluctuations-fluid momentum 
flux coupling induces a gradient in the stress which leads to thermal migration in both systems. 

PACS numbers: 



I. INTRODUCTION 

The first observations of the migration of a species due 
to a temperature gradient were reported by Ludwig and 
Soret more than 100 years ago [3, Since that time, 
the effect has been observed in many multicomponent 
systems including fluid mixtures, colloidal suspensions, 
and polymer solutions (J| . The mass flux of a species is 
described by Fick's law with an added term to account 
for thermal diffusion: 



dc 



dT 



Jy = —pD- pDtc{1 — c) — 



dy 



dy 



(1) 



where Jy is the particle flux in the y-direction. The first 
term denotes diffusion due to a density gradient: D is 
the molecular diffusion coefficient, c is the mass fraction 
of the migrating species, and p is the mass density. The 
second term describes diffusion due to the temperature 
gradient: Dt is the thermal diffusion coefficient and T is 
the temperature. At steady state, Jy = and the Soret 
coefficient, St is defined: 



St = 



Dt 
~D 



1 



dc/dy 



c(l-c) dT/dy' 



(2) 



Note that St can be positive or negative depending on 
whether the species migrates to the hot [St < 0) or cold 
{St > 0) region. 

In general, the thermal diffusion coefficient for a 
molecule will depend on many factors For example, 
a recent simulation conducted on DNA tightly confined 
to nanometer scale channels showed migration towards 
the heated region Q. In contrast, experiments on DNA 
unconfined in the direction of the temperaturegradient 
show the polymer migrating to the cold region 1] . Other 
factors shown in experiment to be important include the 
average temperature of the system, the solvent used, and 
electrostatic effects [1, 0, Hi • 

Since many factors can play an important role in deter- 
mining the thermal diffusion coefficient, theoretical pre- 
dictions and experimental data sometimes do not agree. 



The experiments in used colloidal particles of many 
different materials and showed that Dt will increase with 
increasing particle diameter, a, when some aqueous solu- 
tions are used to suspend the colloids and decrease with 
increasing a for others. These experiments suggest that 
the surface interactions between the particle and the sol- 
vent play important roles in the particles' thermodiffu- 
sion coefficient. Another experiment using carboxyl mod- 
ified polystyrene particles showed that Dt increases with 
increasing a. The authors propose a model based on local 
equilibrium conditions that predicts Dt will only increase 
with increasing gp^. They observe this trend when the 
magnitude of the temperature gradient is much smaller 
than that in the experiments reported in Q . They there- 
fore surmise that the difference between their results and 
those in [§] are due to diflFerences in gradient magnitude 
and nonlinear effects in large gradients. Another model 
based on volume transport theory has proposed that Dt 
for dilute solutions of large molecules depends only on the 
solvent isobaric thermal expansion coefficient assuming 
that the pressure is uniform throughout the fluid [TTI.[T^. 

Several theoretical studies have proposed that the ther- 
mal diffusion of colloids or polymers is a surface driven 
phenomenon. This approach was first adopted by Ruck- 
enstein [l3] . Piazza and Guarino use this model to qual- 
itatively predict the role of electrostatics in the thermal 
diffusion of charged micelles [8] . In general, these studies 
use the hydrodynamic equations to calculate a pressure 
gradient and/or volume force acting locally on the parti- 
cles that is induced by a non-uniform distribution of sol- 
vent molecules or temperature dependent solvent-particle 
interactionji, HI, Qlll, [H, Others also include a 
macroscopic pressure gradient due to the response of the 
solvent alone to the temperature gradient and this model 
nearly quantitatively reproduces the mobility of polymers 
as measured in experiments[14]. Another approach be- 
gins with the kinetic theory of diffusion in a nonuniform 
temperature field to derive expressions for the Soret coef- 
ficient which allow for both positive and negative values 
of S't M- 
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Here, we present results from a lattice Boltzmann sim- 
ulation of A DNA in a microchannel subjected to a ther- 
mal gradient that quantitatively agrees with published 
experimental values for the thermal diffusion coefficient. 
We show that a non- equilibrium stress develops from par- 
ticle fluctuations. This component of the solvent char- 
acteristics causes thermal migration of the species. We 
also investigate the thermal diffusion of small, spheri- 
cal particles and show that Dt decreases with increasing 
diameter, independent of the magnitude of the thermal 
gradient. 



II. LATTICE BOLTZMANN WITH BROWNIAN 
DYNAMICS SIMULATION. 

We use a simulation based on the lattice Boltzmann 
method (LBM) for the fluid coupled with a worm-like 
chain (WLC) inodel with Brownian dynamics (BD) for 
the polymer [l^, [2^ [2l|, The fundamental quantity 
in the LBM is the fluid velocity distribution function, 
ni(r, t), which describes the fraction of fluid particles 
with a discretized velocity, Ci, at each lattice site[23l[2i[. 
To describe the velocities at each node, a 19 discrete ve- 
locity scheme in three dimensions is used. The veloci- 
ties can be represented by: (0,0,0), (±1,0,0), (0, ±1,0), 
(0,0, ±1), (±1,±1,0), (G,±l,±l), and (±1,0, ±1) and 
have magnitudes Ci — |ci| — 0, 1, or \/2. The maximum 
velocity in the simulation is given by the speed of sound: 
Cs = a/ 1/3 Ax/ Ar where Aa; is the lattice spacing and 
At is the time step. The distributions will be Maxwell- 
Boltzmann at equilibrium and can be represented by a 
second order expansion: 

= pa^^ [1 + (ci . u)/cl + uu : {c,c, - c^I)/(2c^)] (3) 

where u is the local velocity. The coefficients a^* are 
found by satisfying the isotropy condition: 

a"Cj_aCi(3C^~,CiS = c'^iSa^S^S + Sa^Sps + SasSpj) (4) 

i 

where a,/?, 7, and 5 represent the x, y, or z axis. The 
equilibrium conditions for the density p, momentum den- 
sity j, and momentum flux density 11: 



j = PU = E Ci • Tlf (6) 



n = p(uu + c2l)=^<^.CiCi. (7) 

i 

must also be satisfied. 



At each time step, the fluid particles will collide with 
their nearest and next nearest neighbors. The velocity 
distributions will evolve according to: 

ni(r + CiAr,t + At) = nj(r, t) + L.,j[nj{r,t) -n^'(r, i)] 

where L is a collision operator for fluid particle collisions 
such that the fluid always relaxes towards the equilibrium 
distribution. In the limit of small Knudson and Mach 
number, this equation has been shown to be equivalent 
to the Navier Stokes equation [25j . 

Collisions are simplified by transforming the from 
velocity space into the hydrodynamic moments, Mq = 
m • n, where Mg is the g*'' moment of the distribution, 
m is the transformation matrix, and n = (no,ni, ...nig). 
The density, momentum density, momentum flux, and 
the kinetic energy fiux constitute the nineteen moments; 
however, kinetic energy fiux moments conserve energy 
and are considered 'ghost' moments. 

The collision operator is chosen to be a diagonal ma- 
trix with elements Tq^,t^^, where Tq is the char- 
acteristic relaxation time of the moment q. The con- 
served moments, density and momentum, are considered 
to have an infinite relaxation time and thus Tq12 3 — ^■ 
The other moments have a single relaxation time, Ts , as 
in the Bhatanagar-Gross-Krook (BGK) model [2Q|. The 
relaxation time is constrained to be smaller than the fiuid 
momentum diffusion time across the system, < pH^ /rj 
where the shear viscosity, 77 is given by: 77 — pc1{ts ~ 0.5) 
for Ts > 0.5 [25]. In our simulations, we have Tg = 1.0. 

A worm-like chain model is adopted for the A DNA 
in the simulation. The model has been parameterized 
to capture molecular dynamics in bulk solution at T = 
298K, and has been used to accurately predict the diffu- 
sivity of 48.5 kbps YOYO-stained DNA in microchannels 

Each molecule is represented by 11 beads connected 
by 10 worm-like springs. The position and velocity of 
the beads are updated using the explicit Euler method. 

u{t + At) = u{t) + {{t)At/m (9) 

x(t + At)^x{t)+u{t)At (10) 

where f(i) is the total force acting on the bead, u{t) is 
the velocity, and x(i) is the position of the bead with 
mass m at time t. The time step for the beads is At. 
The forces acting on the bead include excluded volume 
effects, the elastic force of the springs, hydrodynamic in- 
teractions with the solvent, and the Brownian motion of 
the particles. 

The excluded volume interactions between segments 
are calculated using a Gaussian excluded volume poten- 
tial that leads to self- avoiding walk statistics: 

f'S-^i*.r.«,V3;|j)-p(^*^) (11, 

where ly — af. is the excluded volume interaction param- 
eter, Nks = 19.8 is the number of Kuhn segments per 
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spring, and S'^ — {Nks/Q)^^ is the characteristic size of III. RESULTS 

the coarse grained beads. 

An experimentally determined force-extension relation A. Thermal migration of A DNA 

is used to calculate the elastic force on a bead 30]: 



2<Jk 



[(1- 



Ifj - Til 

NksO-k 



Ifj -nl 



(12) 



0.12 



This was found when measuring the force-extension re- 
lation for the entire chain. The force-extension relation 
is accurate when Nks 3> 1 and we apply this equation to 
single chain segments of Nks = 19.8. 

The fluid exerts a frictional force on the beads, given 
by: 



'/ 



-C(Up - U/) 



(13) 



where Up is the velocity of the bead, u/ is the velocity 
of the fluid at the bead position, C, = Qirrja is the friction 
coefficient, and a is the hydrodynamic radius of the bead 
p^ . The simulation lattice size. As, is chosen to be 0.5 
/XTO. For our model DNA chain, each bead has a hy- 
drodynamic radius of a = 0.077 fim, or 0.154Aa; [28l. [3]|. 
Since the beads' positions are not limited to the lattice 
where the fluid velocity is well defined, u/ at the position 
of the bead is determined by linearly interpolating the 
velocities of the nearest neighbor lattice sites such that 
Uf ~ X^ien T! ""^iUj. The weighting factors Wi are normal- 
ized and U; represents the fluid velocity at site i. The 
momentum transfer to the bead is Aj = —FfAt/Ax^. 
The bead will also transfer this momentum to the fluid. 
The momentum transfer from the bead to a neighbor site 
i with velocity q is Af^ ~ Wipac^Aj ■ Cq [l9[. 

The beads also undergo Brownian motion. The ther- 
mal fluctuations of the beads will be drawn from a Gaus- 
sian distribution with zero mean and a variance that 
varies with the bead height: cry — 2kBT{y)C,At. Here, 



2{T, 



hot 



cold 



\{Y„^axl'^-y)\+T,oid (14) 



where Ymax is the width of the channel, y is the position 
of the bead in the channel, T^ot is the maximum and 
Tcoid is the minimum temperature in the channel. This 
leads to a saw tooth shape of the temperature plotted as 
a function of y. It should be noted that the temperature 
gradient in the system is only implemented here; all other 
forces and fluid properties are independent of location. 

For this work, 50 polymers were simulated in a con- 
tainer of size 20 /im x 20 /im x 20 /xm. Periodic boundary 
conditions were imposed in all directions unless otherwise 
noted. The time step for the fluid is At = 8.8 x 10^'^s, 
and for the polymer At = 1.72 x 10~^s as calculated using 
T = Tcoid- The total simulation time was 1760 seconds. 
Data was recorded once every 17.6 seconds; the flnal 40 
time steps were used to determine Dt- 




FIG. 1: Number fraction, n/ritotai where n is the number of 
beads whose center of mass is between y — 0.25 and y + 0.25 
and ritotai is the total number of beads, as a function of height 
for A-DNA subjected to different temperature gradients. Data 
are averaged over the final 40 time steps of five simulations 
started from different random initial conditions. Shown is the 
average of the two mirror-image halves of the periodic system. 

As shown in Figure [U the DNA accumulates in the 
center of the channel where the temperature minimum is 
found. The simulations with a larger temperature gradi- 
ent result in a larger concentration gradient. The profile 
is nearly flat when the temperature is uniform in the 
channel. These results are in qualitative agreement with 
the work of Duhr and Braun [l| . 

For quantitative comparison, we use the same equation 
found in [l[: 



c(z) 
Co 



exp [~STiT{z) - To)] 



(15) 



where c{z) is the concentration of DNA at position z, cq 
is the maximum concentration, T(z) is the temperature 
at z, To is the temperature at the position of cq, and St is 
the Soret coefficient. This equation is derived from Eqn. 
([T]) with Jy = and c <C 1. The average Soret coefficient 
was found by fitting the density profile to Eqn. (fT51) and 
solving for St- As in [l|, we use D = l.O/im^/s for A 
DNA to calculate Dt from St Hll. We find, for the 
data presented in Figure [U Dt = 0.38±0.1 ^m^/sK. 
This is in good agreement with the value of 0.4 fim? / sK 
reported in [l|. 

In [l| , identical thermal diffusion coefficients were mea- 
sured for 27bp and 48.5 Kbp DNA. We find similar values 
of Dt for 48.5 Kbp {Dt = 0.40±0.06 fim'^/sK), 19.4 
Kbp {Dt = A6±0. 06 ijm^/sK), and 67.9 Kbp DNA 
{Dt = 0.40±0.06 fim^/sK) for AT = iK. The diffusion 
coefficient D of the individual molecules were calculated 
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according to: 

Dl=D,{^)-°-'^' (16) 

where D\ is 1 iJ,m?/s [38], L is the length of the DNA, 
Lx is 48.5 Kbp, the length of A DNA, and Dl is the 
molecular diffusion coefficient for DNA of length L[28j. 
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FIG. 2: Number fraction, n/ritotai wliere n is tlie number of 
beads whose center of mass is between y— 0.25 and J/+0.25 and 
ntotai is the total number of beads, as a function of height at 
different times for A-DNA with AT — 4K. Data are averaged 
over 10 time steps of five simulations started from different 
random initial conditions. Shown is the average of the two 
mirror-image halves of the periodic system. 

The development of the density profile over time can be 
seen in Fig. [2l The profile develops slowly and only af- 
ter more than 300 seconds does the system reach steady 
state. This time frame is similar to that observed by 
Duhr, Arduini, and Braun who report that several hun- 
dred seconds are needed to reach the steady state [l|. 

B. Mechanism of thermal migration 

To understand why the polymers migrate to the colder 
regions, we investigated the dynamics of the solvent. 
Since the properties of the fluid were kept constant across 
the channel, the migration must result from interactions 
between the fluid and the polymers. Thus the momen- 
tum flux of the fluid within 2 lattice sites of a bead was 
recorded. This quantity is known to be coupled to the 
DNA fluctuations [l^ and will contribute to the local 
fluid stress[3^]. A gradient in Hyy = X)i "i^ ' ^iCi is ob- 
served for non-zero temperature gradient but is absent 
for simulations with uniform temperature as seen in Fig- 
ure [31 As gravity is absent, the stress is only due to the 
momentum flux [s^l- This gradient in flux is therefore a 
gradient in stress which causes the polymers to migrate 
into the cold regions. 

The difference in stress is due solely to the interac- 
tion of the beads and the fluid. Without the presence 



of the DNA, the fluid will relax back to equilibrium con- 
ditions and therefore uniform stress across the channel. 
In fact, when the force the polymers exert on the fluid 
is set to zero in the simulation, no thermal diffusion is 
observed. Thus the fluctuations of the polymers induce 
a local, short-lived gradient in stress which causes the 
thermal migration of the species. 

The steady state is reached when the particle flux from 
the cold region to the hot region equals the flux in the 
reverse direction. The particles will migrate only because 
of Brownian motion in this simulation. Thus there is a 
well defined probability of a particle in the cold side re- 
ceiving a sufficient Brownian kick to move to a higher 
stress region. In the steady state, this probability times 
the number of beads in the cold side must equal the cor- 
responding probability a particle will move towards the 
low stress area times the number of beads in the hot re- 
gion. Since the Brownian force depends on the viscosity 
of the fluid and the hydrodynamic radius of the particle, 
the thermal diffusion coefficient should also depend on 
these parameters. We will investigate these dependen- 
cies within the limits of our simulation. 



C. Particle size effects 

The polymer and fluid are coupled through the Brow- 
nian and viscous forces. Both terms depend on ^ = 
67r?7a; the viscous force explicitly and the Brownian term 
through the standard deviation of the distribution of 
the fluctuations. In the simulation, neither the hydro- 
dynamic radius, a, or viscosity of the fluid, rj, appears 
singly. Thus, doubling the value of one term is analogous 
to doubling the other. However, we consider changes to 
C to be changes to a since experimental work has been 
conducted on the effects of changing colloidal particle di- 
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FIG. 3: Stress in the y-direction across the channel for A- 
DNA. Data are averaged over the final 40 time steps of five 
simulations started from different random initial conditions. 
Shown is the average of the two mirror-image halves of the 
periodic system. 
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ameter rather than changing the viscosity of the fluid 

The simulations were conducted using 100 individual 
spheres not attached to each other by springs. We inves- 
tigated particles with a hydrodynamic radius of 0.0385, 
0.077, 0.154, and 0.231 /im with a temperature differ- 
ence of 4K across the lOfim channel. The diffusion co- 
efficient, D, of each species was calculated according to 
D — kBTcoid/^'^'n^ where a is the hydrodynamic radius 
of the particle and rj is the viscosity of the solution. Sim- 
ilar values were obtained for the Soret coefficient, S't, 
as can be seen in Figure IH Here, the density profiles 
are nearly identical for all diameters. However, because 
of the different diffusion coefficients, D, larger values for 
Dt were obtained for the smaller particles. See Table U 
for a comparison of values. It has been suggested that 
Dt decreases with increasing a because the gradient is 
too steep to allow the particles to reach local equilib- 
rium. We therefore decreased the temperature difference 
to 2K to test this hypothesis. However, as Table U shows, 
the values of Dt do not change appreciably. This sug- 
gests the decrease of Dt with increasing radius is not 
due to a particles being farther out of local equilibrium. 
Indeed, all of the above simulations meet the criteria, 
(oS't)^^ > VT, as set by the local equilibrium condition 
in [10]. 

We suggest a different explanation based on the analy- 
sis of the non-equilibrium stress induced by the thermal 
gradient. The fluid stress gradient can be seen in Figure 
O The larger diameter particles induce a steeper gradient 
in the stress than do the smaller spheres. However, the 
steady state density profile is the same for all particles 
due to the differences in the Brownian force on the par- 
ticles. Since the standard deviation of the distribution is 
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FIG. 4: Number fraction, n/ntotai where n is the number of 
beads whose center of mass is between j/— 0.25 and J/-I-0.25 and 
ntotai is the total number of beads, as a function of height for 
spherical particles of different diameter with AT = AK. Data 
are averaged over the final 40 time steps of five simulations 
started from different random initial conditions. Shown is the 
average of the two mirror-image halves of the periodic system. 



TABLE I: Values of the thermal diffusion coefficient, Dt, 
given in ^m? / sK for spheres of different hydrodynamic ra- 
dius, a, and temperature difi'erence, AT across the 10 /im 
channel. 



AT 


a = 0.0385Atm 


a — 0.077/im 


a = 0.154/im 


4K 


2.1 ±0.3 


1.12 ±0.06 


0.59 ±0.04 


2K 


2.3 ±0.4 


1.1 ±0.2 


0.60 ±0.01 



proportional to the hydrodynamic radius, larger particles 
experience larger fluctuations. This is exactly balanced 
by the steeper gradient induced, leading to steady state 
density profiles that are nearly the same for all particle 
sizes. 



IV. CONCLUSIONS. 

We use a lattice-Boltzmann simulation to investigate 
the mechanism behind thermal diffusion of A-DNA and 
spherical particles. This method allows us to capture the 
non-equilibrium stress in the fluid due to the temperature 
gradient. We find good agreement with experimental re- 
sults for the thermal diffusion coefficient, Dt, for A-DNA 
[l[. It is also observed that Dt decreases as the diame- 
ter, a, of diffusing spherical particles increases in partial 
agreement with experiments by Shiunduh, et al. [9|. 

The thermal diffusion coefficient is observed to de- 
crease if the size of the diffusing species is increased. In 
our work, unlike in experiments, the fluid characteris- 
tics such as viscosity are held constant across the chan- 
nel. It is also noted that decreasing the temperature 
difference across the channel did not change these re- 
sults. Therefore the dependence on particle size can not 
be explained as a result of the spheres not reaching lo- 
cal equilibrium with the fluid or the fluid characteristics 
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FIG. 5: Stress in the y-direction across the channel for spheri- 
cal particles of different diameter. Data are averaged over the 
final 40 time steps of five simulations started from different 
random initial conditions. Shown is the average of the two 
mirror-image halves of the periodic system. 
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varying too much across the channel. 

Instead, the non-equiUbrium component of the fluid 
flow is observed to be linked with the migration of the 
solute. The induced stress is significantly more in the 
hot region than in the cold. Thus particles will migrate 
to the cold regions. The thermal diffusion coefficient will 
therefore depend on factors influencing the interaction 
between the solvent and solute. 

This picture is in agreement with the theoretical work 
of several authors [1, [H, [H, [H, [H, [13 who investigate 
local pressure gradients induced by non-isotropic interac- 
tions between the solute and solvent. However, in those 
studies, the specifics of the interaction of the diffusing 
particle and surrounding fluid played an important role 



in determining St- Here, we have quantitatively pre- 
dicted the thermal diffusion coefficient of DNA as well 
as the time scales of the phenomenon without including 
the details of the interaction between the polymer and 
solvent. Instead of these characteristics being important, 
only the Brownian motion of the particle induces a local 
stress gradient in the fluid. This may indicate that hydro- 
dynamic memory effects, already shown to be important 
in the mass diffusion of polymers and colloids, should 
be considered when developing theories of thermal diffu- 
sion [H, [33, . We therefore hope that this work 
will inspire new directions in theoretical examinations of 
thermal diffusion. 
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